home *** CD-ROM | disk | FTP | other *** search
/ Games of Daze / Infomagic - Games of Daze (Summer 1995) (Disc 1 of 2).iso / x2ftp / msdos / math / apm / muldiv.c < prev    next >
Encoding:
Text File  |  1988-10-07  |  15.6 KB  |  564 lines
  1. /******************************************************************************
  2.  
  3.     Arbitrary Precision Math Library General Public License
  4.             (Written October 5, 1988)
  5.  
  6.  Copyright (C) 1988 Lloyd Zusman, Master Byte Software, Los
  7.  Gatos, California.  Everyone is permitted to copy and distribute
  8.  verbatim copies of this license, but changing it is not allowed.
  9.  You can also use this wording to make the terms for other programs.
  10.  
  11.  The wording of this license is based on that of the
  12.  "GNU EMACS GENERAL PUBLIC LICENSE" by Richard Stallman,
  13.  Copyright (C) 1985, 1987, 1988, version of February 11, 1988,
  14.  but since some of the text has been changed, please be sure to
  15.  READ THIS CAREFULLY!
  16.  
  17.   This general public license is intended to give everyone the right
  18. to share the Arbitrary Precision Math Library (hereinafter referred to
  19. as the "APM Library").  To make sure that you get the rights we want
  20. you to have, I need to make restrictions that forbid anyone to deny
  21. you these rights or to ask you to surrender the rights.
  22.  
  23.   Specifically, we want to make sure that you have the right to give
  24. away copies of the APM Library, that you receive source code or else
  25. can get it if you want it, that you can change the APM Library or use
  26. pieces of it in new programs, and that you know you can do these
  27. things.
  28.  
  29.   To make sure that everyone has such rights, we have to forbid you to
  30. deprive anyone else of these rights.  For example, if you distribute
  31. copies of the APM Library, you must give the recipients all the
  32. rights that you have.  You must make sure that they, too, receive or
  33. can get the source code.  And you must tell them their rights.
  34.  
  35.   Also, for our own protection, we must make certain that everyone
  36. finds out that there is no warranty for the APM Library.  If the APM
  37. Library is modified by someone else and passed on, we want its
  38. recipients to know that what they have is not what we distributed, so
  39. that any problems introduced by others will not reflect on our
  40. reputation.
  41.  
  42.   Therefore we (Lloyd Zusman and Master Byte Software) make the
  43. following terms which say what you must do to be allowed to
  44. distribute or change the APM Library.
  45.  
  46.             COPYING POLICIES
  47.  
  48. 1. You may copy and distribute verbatim copies of the APM Library
  49. source code as you receive it, in any medium, provided that you
  50. conspicuously and appropriately publish on each copy a valid copyright
  51. notice "Copyright (C) 1988 Lloyd Zusman, Master Byte Software, Los
  52. Gatos, California" (or with whatever year is appropriate); keep intact
  53. the notices on all files that refer to this License Agreement and to
  54. the absence of any warranty; and give any other recipients of the the
  55. APM Library program a copy of this License Agreement along with the
  56. program.  You may charge a distribution fee for the physical act of
  57. transferring a copy.
  58.  
  59.   2. You may modify your copy or copies of the APM Library source code or
  60. any portion of it, and copy and distribute such modifications under
  61. the terms of Paragraph 1 above, provided that you also do the following:
  62.  
  63.     a) cause the modified files to carry prominent notices stating
  64.     that you changed the files and the date of any change; and
  65.  
  66.     b) cause the whole of any work that you distribute or publish, that in
  67.     whole or in part contains or is a derivative of the APM Library or any
  68.     part thereof, to be licensed to all third parties on terms identical
  69.     to those contained in this License Agreement (except that you may
  70.     choose to grant more extensive warranty protection to some or all
  71.     third parties, at your option).
  72.  
  73.     c) You may charge a distribution fee for the physical act of
  74.     transferring a copy, and you may at your option offer warranty
  75.     protection in exchange for a fee.
  76.  
  77.     d) You may not charge a license fee for the whole of any work that
  78.     you distribute or publish, that in whole or in part contains or is
  79.     a derivative of the APM library or any part thereof, without the
  80.     express written permission of Lloyd Zusman and Master Byte Software;
  81.     whether this permission is granted for free or in return for goods
  82.     services, royalties, or other compensation will be determined
  83.     solely by Lloyd Zusman and Master Byte Software.
  84.  
  85. Mere aggregation of another unrelated program with this program (or its
  86. derivative) on a volume of a storage or distribution medium does not bring
  87. the other program under the scope of these terms.
  88.  
  89.   3. You may copy and distribute the APM Library (or a portion or
  90. derivative of it, under Paragraph 2) in object code or executable form
  91. under all the terms of Paragraphs 1 and 2 above provided that you also
  92. do one of the following:
  93.  
  94.     a) accompany it with the complete corresponding machine-readable
  95.     source code, which must be distributed under the terms of
  96.     Paragraphs 1 and 2 above; or,
  97.  
  98.     b) accompany it with a written offer, valid for at least three
  99.     years, to give any third party free (except for a nominal
  100.     shipping charge) a complete machine-readable copy of the
  101.     corresponding source code, to be distributed under the terms of
  102.     Paragraphs 1 and 2 above; or,
  103.  
  104.     c) accompany it with the information you received as to where the
  105.     corresponding source code may be obtained.  (This alternative is
  106.     allowed only for noncommercial distribution and only if you
  107.     received the program in object code or executable form alone.)
  108.  
  109. For an executable file, complete source code means all the source code
  110. for all modules it contains; but, as a special exception, it need not
  111. include source code for modules which are standard libraries that
  112. accompany the operating system on which the executable file runs.
  113.  
  114.   4. You may not copy, sublicense, distribute or transfer the APM
  115. Library except as expressly provided under this License Agreement.
  116. Any attempt otherwise to copy, sublicense, distribute or transfer the
  117. APM Library is void and your rights to use the APM Library under this
  118. License agreement shall be automatically terminated.  However, parties
  119. who have received computer software programs from you with this
  120. License Agreement will not have their licenses terminated so long as
  121. such parties remain in full compliance.
  122.  
  123.   5. If you wish to incorporate parts of the APM Library into other
  124. programs whose distribution conditions are different, write to Lloyd
  125. Zusman at Master Byte Software.  We have not yet worked out a simple
  126. rule that can be stated here, but we will often permit this.  We will
  127. be guided by the goals of (1) preserving the free status of all
  128. derivatives of our free software; of (2) promoting the sharing and
  129. reuse of software; and of (3) not allowing anyone to profit from the
  130. use of our software without us also having the opportunity to share
  131. in these profits.
  132.  
  133. Your comments and suggestions about our licensing policies and our
  134. software are welcome!  Please contact Lloyd Zusman, Master Byte
  135. Software, 127 Wilder Ave., Los Gatos, California 95030, or call
  136. (408) 395-5693.
  137.  
  138.                NO WARRANTY
  139.  
  140.   BECAUSE THE APM LIBRARY IS LICENSED FREE OF CHARGE, WE PROVIDE
  141. ABSOLUTELY NO WARRANTY, TO THE EXTENT PERMITTED BY APPLICABLE STATE
  142. LAW.  EXCEPT WHEN OTHERWISE STATED IN WRITING, MASTER BYTE SOFTWARE,
  143. LLOYD ZUSMAN AND/OR OTHER PARTIES PROVIDE THE APM LIBRARY "AS IS"
  144. WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING,
  145. BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
  146. FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY
  147. AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE THE APM
  148. LIBRARY PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY
  149. SERVICING, REPAIR OR CORRECTION.
  150.  
  151.   IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW WILL MASTER BYTE
  152. SOFTWARE, LLOYD ZUSMAN, AND/OR ANY OTHER PARTY WHO MAY MODIFY AND
  153. REDISTRIBUTE THE APM LIBRARY AS PERMITTED ABOVE, BE LIABLE TO YOU FOR
  154. DAMAGES, INCLUDING ANY LOST PROFITS, LOST MONIES, OR OTHER SPECIAL,
  155. INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR
  156. INABILITY TO USE (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA
  157. BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY THIRD PARTIES OR A
  158. FAILURE OF THE PROGRAM TO OPERATE WITH PROGRAMS NOT DISTRIBUTED BY
  159. MASTER BYTE SOFTWARE) THE PROGRAM, EVEN IF YOU HAVE BEEN ADVISED OF
  160. THE POSSIBILITY OF SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY.
  161.  
  162. ******************************************************************************/
  163.  
  164.  
  165. /*
  166.  * Multiplication and division routines for the APM library.
  167.  *
  168.  * $Log:    muldiv.c,v $
  169.  * Revision 1.0  88/10/05  12:38:15  ljz
  170.  * Initial release.
  171.  * 
  172.  */
  173. #ifndef lint
  174. static char rcsid[] = "$Header: muldiv.c,v 1.0 88/10/05 12:38:15 ljz Exp $";
  175. #endif /* ! lint */
  176.  
  177. #include <stdio.h>
  178. #include "apm.h"
  179. #include "apmlocal.h"
  180.  
  181. int
  182. apm_multiply(result, num1, num2)
  183. APM result;
  184. APM num1;
  185. APM num2;
  186. {
  187.     int length1;
  188.     int length2;
  189.     int len;
  190.     int dp;
  191.     int n1;
  192.     int n2;
  193.     int sum;
  194.     short sign;
  195.     short base;
  196.     short carry;
  197.     long tempval;
  198.  
  199.     apm_errno = APM_OK;
  200.  
  201.     ERR_RETURN(APM_val_format(result));
  202.     ERR_RETURN(apm_validate(num1));
  203.     ERR_RETURN(apm_validate(num2));
  204.  
  205.     if (num1->base != num2->base) {
  206.         return (APM_error(APM_EBASE));
  207.     }
  208.  
  209.     if (result == num1 || result == num2) {
  210.         return (APM_error(APM_EOVERLAP));
  211.     }
  212.  
  213.     base = num1->base;
  214.  
  215.     sign = SIGNOF(num1->sign) * SIGNOF(num2->sign);
  216.  
  217.     ERR_RETURN(APM_trim(num1, 1, 1));
  218.     ERR_RETURN(APM_trim(num2, 1, 1));
  219.  
  220.     length1 = num1->length;
  221.     length2 = num2->length;
  222.     
  223.     len = length1 + length2;
  224.     dp = num1->dp + num2->dp;
  225.  
  226.     ERR_RETURN(APM_size(result, len));
  227.  
  228.     APM_zero_shorts(result->data, len);
  229.  
  230.     for (n1 = 0; n1 < length1; ++n1) {
  231.         carry = 0;
  232.         sum = n1;
  233.         for (n2 = 0; n2 < length2; ++n2, ++sum) {
  234.             tempval = num1->data[n1];
  235.             tempval *= num2->data[n2];
  236.             tempval += result->data[sum] + carry;
  237.             result->data[sum] =
  238.                 (short)(tempval % base);
  239.             carry =
  240.                 (short)(tempval / base);
  241.         }
  242.         result->data[sum] = carry;
  243.     }
  244.  
  245.     result->length = len;
  246.     result->dp = dp;
  247.     result->base = base;
  248.     result->sign = sign;
  249.  
  250.     return (APM_error(APM_trim(result, 1, 1)));
  251. }
  252.  
  253. /*
  254.  * For division, I use the algorithm described in Knuth, "The Art of
  255.  * Computer Programming, Volume 2, Seminumerical Algorithms", second
  256.  * edition, pp. 255-260.
  257.  */
  258. int
  259. apm_divide(result, frac_precision, remainder, num1, num2)
  260. APM result;
  261. int frac_precision;
  262. APM remainder;
  263. APM num1;
  264. APM num2;
  265. {
  266.     static APM temp = (APM)NULL;
  267.     static short *dividend = (short *)NULL;
  268.     static int dividendLength = 0;
  269.     static short *divisor = (short *)NULL;
  270.     static short *xdivisor = (short *)NULL;
  271.     static int divisorLength = 0;
  272.  
  273.     int i;
  274.     int j;
  275.     int n;
  276.     int uJ;
  277.     int vJ;
  278.     int M;
  279.     int N;
  280.     int MplusN;
  281.     int Nplus1;
  282.     int MplusNplus1;
  283.     int length1;
  284.     int length2;
  285.     int offset;
  286.     int left;
  287.     int frac;
  288.     short divN1;
  289.     short divN2;
  290.     short sign;
  291.     short base;
  292.     short scaleFactor;
  293.     short borrow;
  294.     short qHat;
  295.     long temp1;
  296.     long temp2;
  297.     long temp3;
  298.  
  299.     apm_errno = APM_OK;
  300.  
  301.     if (frac_precision <= 0) {
  302.         return (APM_set_errno(APM_EPARM));
  303.     }
  304.  
  305.     ERR_RETURN(APM_val_format(result));
  306.     if (remainder != NULL) {
  307.         ERR_RETURN(APM_val_format(remainder));
  308.     }
  309.     ERR_RETURN(apm_validate(num1));
  310.     ERR_RETURN(apm_validate(num2));
  311.  
  312.     if (result == remainder || result == num1 || result == num2) {
  313.         return (APM_error(APM_EOVERLAP));
  314.     }
  315.     if (remainder == num1 || remainder == num2) {
  316.         return (APM_error(APM_EOVERLAP));
  317.     }
  318.  
  319.     base = num1->base;
  320.  
  321.     ERR_RETURN(APM_trim(num2, 1, 1));
  322.     if (num2->length <= 0) {
  323.         /*
  324.          * Division by zero.
  325.          */
  326.         result->length = 0;
  327.         result->dp = 0;
  328.         result->base = base;
  329.         if (remainder != (APM)NULL) {
  330.             remainder->length = 0;
  331.             remainder->dp = 0;
  332.             remainder->base = base;
  333.         }
  334.         return (APM_error(APM_WDIVBYZERO));
  335.     }
  336.  
  337.     ERR_RETURN(APM_trim(num1, 1, 1));
  338.     if (num1->length <= 0) {
  339.         /*
  340.          * Dividend is zero, so result is zero.
  341.          */
  342.         result->length = 0;
  343.         result->dp = 0;
  344.         result->base = base;
  345.         if (remainder != (APM)NULL) {
  346.             remainder->length = 0;
  347.             remainder->dp = 0;
  348.             remainder->base = base;
  349.         }
  350.         return (APM_OK);
  351.     }
  352.  
  353.     sign = num1->sign * num2->sign;
  354.  
  355.     length1 = num1->length;
  356.     length2 = num2->length;
  357.     for (offset = length2; offset > 0; --offset) {
  358.         if (num2->data[offset - 1] != 0) {
  359.             break;
  360.         }
  361.     }
  362.  
  363.     offset = length2 - offset;
  364.  
  365.     frac = frac_precision;
  366.     if (base == SPECIAL_BASE) {
  367.         frac = (frac + (SPECIAL_SCALE - 1)) / SPECIAL_SCALE;
  368.     }
  369.  
  370.     
  371.     left = (length1 - num1->dp) - (length2 - num1->dp);
  372.     if (left < 0) {
  373.         left = 0;
  374.     }
  375.  
  376.     M = frac + left + 1;
  377.     if (M < length1) {
  378.         M = length1;
  379.     }
  380.  
  381.     N = (length2 - offset) + 1;
  382.     Nplus1 = N + 1;
  383.     MplusN = M + N;
  384.     MplusNplus1 = MplusN + 1;
  385.  
  386.     if (temp == (APM)NULL) {
  387.         temp = APM_alloc();
  388.         if (temp == (APM)NULL) {
  389.             return (APM_error(APM_ENOMEM));
  390.         }
  391.     }
  392.  
  393.     ERR_RETURN(APM_size(temp, M + 1));
  394.  
  395.     if (MplusNplus1 > dividendLength) {
  396.         int xlen = MplusNplus1;
  397.         if (xlen < 8) {
  398.             xlen = 8;
  399.         }
  400.         if (dividendLength < 1) {
  401.             dividend = (short *)APM_alloc_mem(NULL, xlen,
  402.                         sizeof (short));
  403.         }
  404.         else {
  405.             dividend = (short *)APM_alloc_mem(dividend, xlen,
  406.                         sizeof (short));
  407.         }
  408.         if (dividend == (short *)NULL) {
  409.             return (APM_error(APM_ENOMEM));
  410.         }
  411.         dividendLength = xlen;
  412.     }
  413.  
  414.     if (Nplus1 > divisorLength) {
  415.         int xlen = Nplus1;
  416.         if (xlen < 8) {
  417.             xlen = 8;
  418.         }
  419.         if (divisorLength < 1) {
  420.             divisor = (short *)APM_alloc_mem(NULL, xlen,
  421.                               sizeof (short));
  422.         }
  423.         else {
  424.             divisor = (short *)APM_alloc_mem(divisor, xlen,
  425.                               sizeof (short));
  426.         }
  427.         if (divisor == (short *)NULL) {
  428.             return (APM_error(APM_ENOMEM));
  429.         }
  430.         if (divisorLength < 1) {
  431.             xdivisor = (short *)APM_alloc_mem(NULL, xlen,
  432.                                sizeof (short));
  433.         }
  434.         else {
  435.             xdivisor = (short *)APM_alloc_mem(xdivisor, xlen,
  436.                                sizeof (short));
  437.         }
  438.         if (xdivisor == (short *)NULL) {
  439.             return (APM_error(APM_ENOMEM));
  440.         }
  441.         divisorLength = xlen;
  442.     }
  443.  
  444.     i = MplusN - length1;
  445.  
  446.     APM_zero_shorts(dividend, i);
  447.     APM_copy_shorts(÷nd[i], num1->data, length1);
  448.  
  449.     divisor[0] = 0;
  450.     APM_copy_shorts(&divisor[1], num2->data, length2 - offset);
  451.  
  452.     /*
  453.      * Knuth: step D1.
  454.      */
  455.     scaleFactor = base / (divisor[N - 1] + 1);
  456.     if (scaleFactor <= 1) {
  457.         dividend[MplusN] = 0;
  458.     }
  459.     else {
  460.         dividend[MplusN] = APM_scalar_mul(dividend, dividend,
  461.                                MplusN,
  462.                                scaleFactor, base);
  463.         APM_scalar_mul(divisor, divisor, N, scaleFactor, base);
  464.     }
  465.  
  466.     divN1 = divisor[N - 1];
  467.     divN2 = divisor[N - 2];
  468.  
  469.     /*
  470.      * Knuth: steps D2 and D7.
  471.      */
  472.     for (j = 0, vJ = M, uJ = MplusN; j <= M; ++j, --vJ, --uJ) {
  473.         /*
  474.          * Knuth: step D3.
  475.          */
  476.         temp1 = dividend[uJ];
  477.         temp1 *= base;
  478.         temp1 += dividend[uJ - 1];
  479.  
  480.         if (dividend[uJ] == divN1) {
  481.             qHat = base - 1;
  482.         }
  483.         else {
  484.             temp2 = temp1 / divN1;
  485.             qHat = (short)temp2;
  486.         }
  487.  
  488.         for (;;) {
  489.             temp2 = divN2;
  490.             temp2 *= qHat;
  491.  
  492.             temp3 = divN1;
  493.             temp3 *= qHat;
  494.             temp3 = temp1 - temp3;
  495.             temp3 *= base;
  496.             temp3 += dividend[uJ - 2];
  497.  
  498.             if (temp2 > temp3) {
  499.                 --qHat;
  500.             }
  501.             else {
  502.                 break;
  503.             }
  504.         }
  505.  
  506.         /*
  507.          * Knuth: step D4.
  508.          */
  509.         xdivisor[N] = APM_scalar_mul(xdivisor, divisor, N, qHat,
  510.                           base);
  511.         borrow = APM_array_sub(&(dividend[vJ]), xdivisor,
  512.                    Nplus1, base);
  513.  
  514.         /*
  515.          * Knuth: step D5.
  516.          */
  517.         if (borrow == 0) {
  518.             temp->data[vJ] = qHat;
  519.         }
  520.         else {
  521.             /*
  522.              * Knuth: step D6.
  523.              */
  524.             temp->data[vJ] = qHat - 1;
  525.             divisor[N] = 0;
  526.             APM_array_add(&(dividend[vJ]), divisor, Nplus1, base);
  527.         }
  528.     }
  529.  
  530.     /*
  531.      * Knuth: step D8.
  532.      */
  533.     temp->length = M + 1;
  534.     temp->dp = (M - offset) -
  535.         ((length1 - num1->dp) - (length2 - num2->dp));
  536.     temp->base = base;
  537.     temp->sign = sign;
  538.  
  539.     ERR_RETURN(apm_round(result, temp, frac_precision));
  540.     ERR_RETURN(APM_trim(result, 1, 1));
  541.  
  542.     if (remainder != (APM)NULL) {
  543.         ERR_RETURN(APM_size(temp, length2));
  544.         if (scaleFactor > 1) {
  545.             ERR_RETURN(APM_scalar_div(temp, ÷nd[1],
  546.                              length2,
  547.                              num2->dp,
  548.                              num2->sign,
  549.                              base, scaleFactor));
  550.         }
  551.         else {
  552.             temp->length = length2;
  553.             temp->dp = num2->dp;
  554.             temp->sign = num2->sign;
  555.             temp->base = base;
  556.             APM_copy_shorts(temp->data, ÷nd[1], length2);
  557.         }
  558.         ERR_RETURN(apm_round(remainder, temp, frac_precision));
  559.         ERR_RETURN(APM_trim(remainder, 1, 1));
  560.     }
  561.  
  562.     return (APM_OK);
  563. }
  564.